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I Abstract: We study graph estimation and density estimation in high dimensions, 

■ using a family of density estimators based on forest structured undirected graphical 
^ , models. For density estimation, we do not assume the true distribution corresponds 

' to a forest; rather, we form kernel density estimates of the bivariate and univariate 

■ marginals, and apply Kruskal's algorithm to estimate the optimal forest on held out 
\ data. We prove an oracle inequality on the excess risk of the resulting estimator relative 
' to the risk of the best forest. For graph estimation, we consider the problem of estimat- 

I— 1 1 ing forests with restricted tree sizes. We prove that finding a maximum weight spanning 

forest with restricted tree size is NP-hard, and develop an approximation algorithm for 
' this problem. Viewing the tree size as a complexity parameter, we then select a forest 

^ . using data splitting, and prove bounds on excess risk and structure selection consis- 

^ \ tency of the procedure. Experiments with simulated data and microarray data indicate 

I "^i ' that the methods are a practical alternative to Gaussian graphical models. 

\ Keywords and phrases: kernel density estimation, tree structured Markov network, 

^ ' high dimensional inference, risk consistency, structure selection consistency. 
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1. Introduction 



One way to explore the structure of a high dimensional distribution P for a random vector 
X = (Xi, . . . ,Xd) is to estimate its undirected graph. The undirected graph G associated 
with P has d vertices corresponding to the variables Xi, . . . , X^, and omits an edge between 
two nodes Xi and Xj if and only if X^ and Xj are conditionally independent given the 
other variables. Currently, the most popular methods for estimating G assume that the 
distribution P is Gaussian. Finding the graphical structure in this case amounts to estimating 
the inverse covariance matrix Q; the edge between Xj and X^ is missing if and only if Qji^ = 0. 
Algorithms for optimizing the £i-regularized log-likelihood have recently been proposed that 
efficiently produce sparse estimates of the inverse covariance matrix and the underlying graph 
(Banerjee, El Ghaoui and d'Aspremont, 2008; Friedman, Hastie and Tibshirani, 2007). 

In this paper our goal is to relax the Gaussian assumption and to develop nonparametric 
methods for estimating the graph of a distribution. Of course, estimating a high dimensional 
distribution is impossible without making any assumptions. The approach we take here is 
to force the graphical structure to be a forest, where each pair of vertices is connected by 
at most one path. Thus, we relax the distributional assumption of normality but we restrict 
the family of undirected graphs that are allowed. 

If the graph for P is a forest, then a simple conditioning argument shows that its density 
p can be written as 

-J 
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where E is the set of edges in the forest (Lauritzen, 1996). Here p{xi,Xj) is the bivariate 
marginal density of variables Xi and Xj, and p{xk) is the univariate marginal density of 
the variable X^. With this factorization, we see that it is only necessary to estimate the 
bivariate and univariate marginals. Given any distribution P with density p, there is a tree 
T and a density pr whose graph is T and which is closest in KuUback-Leibler divergence to 
p. When P is known, then the best fitting tree distribution can be obtained by Kruskal's 
algorithm (Kruskal, 1956), or other algorithms for finding a maximum weight spanning tree. 
In the discrete case, the algorithm can be applied to the estimated probability mass func- 
tion, resulting in a procedure originally proposed by Chow and Liu (1968). Here we are 
concerned with continuous random variables, and we estimate the bivariate marginals with 
nonparametric kernel density estimators. 

In high dimensions, fitting a fully connected spanning tree can be expected to overfit. We 
regulate the complexity of the forest by selecting the edges to include using a data splitting 
scheme, a simple form of cross validation. In particular, we consider the family of forest 
structured densities that use the marginal kernel density estimates constructed on the first 
partition of the data, and estimate the risk of the resulting densities over a second, held 
out partition. The optimal forest in terms of the held out risk is then obtained by finding a 
maximum weight spanning forest for an appropriate set of edge weights. 

A closely related approach is proposed by Bach and Jordan (2003), where a tree is esti- 
mated for the random vector Y = WX instead of X, where is a linear transformation, 
using an algorithm that alternates between estimating W and estimating the tree T. Kernel 
density estimators are used, and a regularization term that is a function of the number of 
edges in the tree is included to bias the optimization toward smaller trees. We omit the 
transformation W, and we use a data splitting method rather than penalization to choose 
the complexity of the forest. 

While tree and forest structured density estimation has been long recognized as a useful 
tool, there has been little theoretical analysis of the statistical properties of such density 
estimators. The main contribution of this paper is an analysis of the asymptotic properties of 
forest density estimation in high dimensions. We allow both the sample size n and dimension 
d to increase, and prove oracle results on the risk of the method. In particular, we assume 
that the univariate and bivariate marginal densities lie in a Holder class with exponent /3 
(see Section 4 for details), and show that 

/ d \\ 

^^^^i^d) L/3/(2+2/3) + ^/3/(l+2/3) I I (l'^) 

where R denotes the risk, the expected negative log-likelihood, k is the number of edges in 
the estimated forest F, and k* is the number of edges in the optimal forest F* that can be 
constructed in terms of the kernel density estimates p. 
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In addition to the above results on risk consistency, we establish conditions under which 

= i^fVl (1-3) 

as n — > oo, where F^^''^ is the oracle forest — the best forest with k edges; this result allows 
the dimensionality d to increase as fast as a (exp(n^/'^^+^))) , while still having consistency in 
the selection of the oracle forest. 

Among the only other previous work analyzing tree structured graphical models is Tan et al. 
(2009) and Chechetka and Guestrin (2007). Tan et al. (2009) analyze the error exponent in 
the rate of decay of the error probability for estimating the tree, in the fixed dimension set- 
ting, and Chechetka and Guestrin (2007) give a PAG analysis. An extension to the Gaussian 
case is given by Tan, Anandkumar and Willsky (2009). 

We also study the problem of estimating forests with restricted tree sizes. In many ap- 
plications, one is interested in obtaining a graphical representation of a high dimensional 
distribution to aid in interpretation. For instance, a biologist studying gene interaction net- 
works might be interested in a visualization that groups together genes in small sets. Such 
a clustering approach through density estimation is problematic if the graph is allowed to 
have cycles, as this can require marginal densities to be estimated with many interacting 
variables. Restricting the graph to be a forest circumvents the curse of dimensionality by 
requiring only univariate and bivariate marginal densities. The problem of clustering the 
variables into small interacting sets, each supported by a tree-structured density, becomes 
the problem of estimating a maximum weight spanning forest with a restriction on the size 
of each component tree. As we demonstrate, estimating restricted tree size forests can also 
be useful in model selection for the purpose of risk minimization. Limiting the tree size gives 
another way of regulating tree complexity that provides larger family of forest to select from 
in the data splitting procedure. 

While the problem of finding a maximum weight forest with restricted tree size may be 
natural, it appears not to have been studied previously. We prove that the problem is NP- 
hard through a reduction from the problem of Exact 3-Gover (Garey and Johnson, 1979), 
where we are given a set X and a family S of 3-element subsets of X, and must choose a 
subfamily of disjoint 3-element subsets to cover X. While finding the exact optimum is hard, 
we give a practical 4-approximation algorithm for finding the optimal tree restricted forest; 
that is, our algorithm outputs a forest whose weight is guaranteed to be at least \w{F*), 
where w{F*) is the weight of the optimal forest. This approximation guarantee translates into 
excess risk bounds on the constructed forest using our previous analysis. Our experimental 
results with this approximation algorithm show that it can be effective in practice for forest 
density estimation. 

In Section 2 we review some background and notation. In Section 3 we present a two-stage 
algorithm for estimating high dimensional densities supported by forests, and we provide a 
theoretical analysis of the algorithm in Section 4, with the detailed proofs collected in an 
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appendix. In Section 5, we explain how to estimate maximum weight forests with restricted 
tree size. In Section 6 we present experiments with both simulated data and gene microarray 
datasets, where the problem is to estimate the gene-gene association graphs. 



2. Preliminaries and Notation 

Let p*{x) be a probability density with respect to Lebesgue measure /i(-) on M.'^ and let 
X^^\ . . . , X*^"^ be n independent identically distributed M''- valued data vectors sampled from 
p*{x) where X^'^ = (xf \ . . . , X^'^). Let Xj denote the range of xf and let = Afi x ■ ■ ■ x X^. 

A graph is a forest if it is acyclic. If F is a d-node undirected forest with vertex set 
Vp = {1, . . . ,d} and edge set E{F) C {1, . . . ,d} x {1, . . . , d}, the number of edges satisfies 
\E{F)\ < d, noting that we do not restrict the graph to be connected. We say that a 
probability density function p{x) is supported by a forest F if the density can be written as 

where each p{xi bivariate density on Xi x Xj, and each p{xk) is a univariate density 

on Xk- More details can be found in Lauritzen (1996). 

Let J-rf be the family of forests with d nodes, and let Vd be the corresponding family of 
densities: 

"^(i = |p > : J p{x) djj,{x) = 1, and p{x) satisfies (2.1) for some F G J^d^ ■ (2.2) 

To bound the number of labeled spanning forests on d nodes, note that each such forest can 
be obtained by forming a labeled tree on c? + 1 nodes, and then removing node d+ 1. From 
Cayley's formula (Aigner and Ziegler, 1998; Cayley, 1889), we then obtain the following. 

Proposition 2.1. The size of the collection J^d of labeled forests on d nodes satisfies 

\J^d\ < {d+lY~\ (2.3) 

Define the oracle forest density 

q* = argmin D(p*|| q) (2.4) 
where the Kullback-Leibler divergence D{p\\q) between two densities p and q is 

D{p\\q)= [ p{x)\og^dx, (2.5) 

Jx q{x) 

under the convention that 01og(0/g) = 0, and p\og{p/0) = oo for p 7^ 0. The following is 
proved by Bach and Jordan (2003). 
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Proposition 2.2. Let q* be defined as in (2.4). There exists a forest F* G Td, such that 

=pf*= n ^%#T n p*(^^) (2-6) 

where p*{xi,Xj) andp*{xi) are the bivariate and univariate marginal densities of p* . 
For any density q{x), the negative log-likeliliood risk R{q) is defined as 

R{q) = -E\ogq{X) = - p* (x) log q{x) dx (2.7) 

Jx 

where the expectation is defined with respect to the distribution of X. 

It is straightforward to see that the density q* defined in (2.4) also minimizes the negative 
log-likehhood loss: 

q* = argmini5(p*|| q) = argmini?(g). (2.8) 
Let p{x) be the kernel density estimate, we also define 

R{q) = — / p{x)\ogq{x) dx. (2.9) 
Jx 

We thus define the oracle risk as R* = R{q*)- Using Proposition 2.2 and equation (2.1), 
we have 

R* = R(q*) = R(p*^^) 

^ ~ zl I P*{xi,Xj)\og ''J dx^dxj- V / p*{xk)\ogp*{xk)dxk 
= - Yl HXf,X,)+ H{Xk), (2.10) 



F* 



where 

HX,, X,) = / fix., X,) log f^i^;f/j dxdx, (2.11) 
is the mutual information between the pair of variables X^, Xj and 

H{Xk) = -! p*{xk)\ogp*{xk)dxk (2.12) 

is the entropy. While the best forest will in fact be a spanning tree, the densities p*{xi,Xj) 
are in practice not known. We estimate the marginals using finite data, in terms of a kernel 
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density estimates Pnj{xi,Xj) over a training set of size rii. With these estimated marginals, 
we consider all forest density estimates of the form 

pf{x)= n -%T^^ n^-(^'^)- (2.13) 

Within this family, the best density estimate may not be supported on a full spanning tree, 
since a full tree will in general be subject to overfitting. Analogously, in high dimensional 
linear regression, the optimal regression model will generally be a full (i- dimensional fit, with 
a nonzero parameter for each variable. However, when estimated on finite data the variance 
of a full model will dominate the squared bias, resulting in overfitting. In our setting of 
density estimation we will regulate the complexity of the forest by cross validating over a 
held out set. 

There are several different ways to judge the quality of a forest structured density esti- 
mator. In this paper we concern ourselves with prediction and structure estimation. 

Definition 2.1 (Risk consistency). For an estimator g„ G Vd, the excess risk is defined 
as R{qn) — R* ■ The estimator g„ is risk consistent with convergence rate 5„ if 

lim limsup P (i?(g„) - R* > M5„) = 0. (2.14) 
In this case we write R{qn) — R* = Op(5„). 

Definition 2.2 (Estimation consistency). An estimator g„ G Vd is estimation consistent 
with convergence rate 6n, with respect to the Kullback-Leibler divergence, if 

lim limsup P(D(p^*|| g„) > M5„) = 0. (2.15) 

Definition 2.3 (Structure selection consistency). An estimator g„ G Vd supported by 
a forest F„ is structure selection consistent if 

P (^E{Fn) + E{F*)^ ^ 0, (2.16) 

as n goes to infinity, where F* is defined in (2.6). 

Later we will show that estimation consistency is almost equivalent to risk consistency. 
If the true density is given, these two criteria are exactly the same; otherwise, estimation 
consistency requires stronger conditions than risk consistency. 

It is important to note that risk consistency is an oracle property, in the sense that the 
true density p*{x) is not restricted to be supported by a forest; rather, the property assesses 
how well a given estimator approximates the best forest density (the oracle) within a class. 
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3. Kernel Density Estimation For Forests 



If the true density p*{x) were known, by Proposition 2.2, tlie density estimation problem 
would be reduced to finding the best forest structure F^, satisfying 

= argmini?(p^) = argmin _D(p*|| p^). (3.1) 

The optimal forest can be found by minimizing the right hand side of (2.10). Since the 
entropy term H{X) = Ylik-^i-^k) is constant across all forests, this can be recast as the 
problem of finding the maximum weight spanning forest for a weighted graph, where the 
weight of the edge connecting nodes i and j is I{Xi]Xj). Kruskal's algorithm (Kruskal, 
1956) is a greedy algorithm that is guaranteed to find a maximum weight spanning tree 
of a weighted graph. In the setting of density estimation, this procedure was proposed by 
Chow and Liu (1968) as a way of constructing a tree approximation to a distribution. At each 
stage the algorithm adds an edge connecting that pair of variables with maximum mutual 
information among all pairs not yet visited by the algorithm, if doing so does not form a 
cycle. When stopped early, after k < d — 1 edges have been added, it yields the best fc-edge 
weighted forest. 

Of course, the above procedure is not practical since the true density p*{x) is unknown. 
We replace the population mutual information I{Xi;Xj) in (2.10) by the plug-in estimate 
In{Xi, Xj), defined as 

Tn{Xi,Xj)= [ Pn{xi,Xj) log ^^"^^'1^/^ dxidxj (3.2) 
where pn{xi,Xj) and Pn{xi) are bivariate and univariate kernel density estimates. Given 



this estimated mutual information matrix M„ 



In{Xi, Xj) 



, we can then apply Kruskal's 



algorithm (equivalently, the Chow-Liu algorithm) to find the best forest structure 

Since the number of edges of Fn controls the number of degrees of freedom in the final 
density estimator, we need an automatic data- dependent way to choose it. We adopt the 
following two-stage procedure. First, randomly partition the data into two sets T>i and 
of sizes ni and ?t,2; then, apply the following steps: 

1. Using "Di, construct kernel density estimates of the univariate and bivariate marginals 
and calculate /„j(Xj,Xj) for i,j G {1, . . . with i ^ j. Construct a full tree F^^ 
with d — 1 edges, using the Chow-Liu algorithm. 

2. Using 1)2, prune the tree Fni'^"* to find a forest Fni with k edges, for < A; < (i — 1. 

Once Fni is obtained in Step 2, we can calculate P'~-(k) according to (2.1), using the kernel 
density estimates constructed in Step 1. 
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3.1. Step 1: Estimating the marginals 



Step 1 is carried out on the dataset Pi. Let K{-) be a univariate kernel function. Given 
an evaluation point the bivariate kernel density estimate for based on the 

observations {x\^\x^^^}s&Vi is defined as 

,,, .^,.1 ^ ^ (^r_^ . ,3,, 

where we use a product kernel with /i2 > be the bandwidth parameter. The univariate 
kernel density estimate PnAxk) for is 




Pn^x,) = -l^Y,^\ \ I , (3.4) 



where /ii > is the univariate bandwidth. Detailed specifications for K{-) and /ii, h2 will be 
discussed in the next section. 

We assume that the data lie in a (i-dimensional unit cube X = [0, 1]*^. To calculate the 
empirical mutual information J„j (Xj, Xj), we need to numerically evaluate a two-dimensional 
integral. To do so, we calculate the kernel density estimates on a grid of points. We choose m 
evaluation points on each dimension, xu < X2i < ■ ■ ■ < Xmi for the ith variable. The mutual 
information J„^(Xj,Xj) is then approximated as 

Tn,{Xi,Xj) = V V {xki , xej ) log ^ -^"^ ^^'Z . (3.5) 

^ fr'^j^^ PnAXki)PnAx£j) 

The approximation error can be made arbitrarily small by choosing m sufficiently large. 
As a practical concern, care needs to be taken that the factors Pmixki) and Pmixej) in the 
denominator are not too small; a truncation procedure can be used to ensure this. Once 
the d X d mutual information matrix M„j = /„^(Xj,Xj) is obtained, we can apply the 
Chow-Liu (Kruskal) algorithm to find a maximum weight spanning tree. 



3.2. Step 2: Optimizing the forest 

The full tree obtained in Step 1 might have high variance when the dimension d is 

large, leading to overfitting in the density estimate. In order to reduce the variance, we prune 
the tree; that is, we choose forest with k < d — 1 edges. The number of edges /c is a tuning 
parameter that induces a bias- variance tradeoff. 

In order to choose k, note that in stage k of the Chow-Liu algorithm we have an edge set 
E^^^ (in the notation of the Algorithm 3.1) which corresponds to a forest F^f^ with k edges, 
where F^i^ is the union of d disconnected nodes. To select k, we choose among the d trees 
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Algorithm 3.1 Chow-Liu (Kruskal) 



1: Input dataPi = 

2: Calculate M„j, according to (3.3), (3.4), and (3.5). 

3: Initialize E^"^ = 

4: for fc = 1 do 

5: (i^'^\i''^'') ■'r- arginax^j J-, M„j(i,j) such that E^^~^'^ U {(i'-'^-', j'-'^-')} does not contain a cycle 

6: £;W ^ u{(i('=),j('=))} 

7: Output tree F^f"^^ with edge set E^'^-^\ 



Let Pn2{xi, Xj) and Pn2(2^fc) be defined as in (3.3) and (3.4), but now evaluated solely based 
on the held-out data in V2. For a density pp that is supported by a forest F, we define the 
held-out negative log-likelihood risk as 

RriiiPp) (3.6) 

^ ~ I Pn2{xi,Xj)\og / dxidxj - y2 / p„2(a;fc)logp(xfc)rfxfc. 

{i,j)eEp p(x, j ^^^^ J^^ 

The selected forest is then F^^'' where 

A;= argmin R^^ (pec^)) (3.7) 
fce{o,...,d-i} ^ "1 ^ 

and where pp(fc) is computed using the density estimate Pni constructed on Vi. 
For computational simplicity, we can also estimate k as 



k = argmax — log 
ke{o,...,d-i} ^gp^ 



\ {i,j)(^EW Pni l-^i ) Pni y^j ) fceV (fc) 



(3.^ 



= argmax — > log , , . (3.9) 

keio,...4-i} n2 VJi.)Pn,(xf))p„,(xf)^ 

This minimization can be efficiently carried out by iterating over the d — 1 edges in Fn'l~^\ 
Once k is obtained, the final forest density estimate is given by 

Pn(x)= n -%%^n^-(^^)- (3.10) 

(i,j)<^EW *= 

Remark 3.1. For computational efficiency. Step 1 can be carried out simultaneously with 
Step 2. In particular, during the Chow-Liu iteration, whenever an edge is added to E^''\ 
the log-likelihood of the resulting density estimator on V2 can be immediately computed. A 
more efficient algorithm to speed up the computation of the mutual information matrix is 
discussed in Appendix B. 
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3.3. Building a forest on held- out data 

Another approach to estimating the forest structure is to estimate the marginal densities on 
the training set, but only build graphs on the held-out data. To do so, we first estimate the 
univariate and bivariate kernel density estimates using Pi, denoted by Pm {xi) and {xi, xj). 
We also construct a new set of univariate and bivariate kernel density estimates using 1)2, 
p„2(xi) and Pn2(^i) Xj)- We then estimate the "cross-entropies" of the kernel density estimates 
for each pair of variables by computing 

In2,ni{Xi,Xj) = / Pn2 i^i, Xj) log ^ ^"^ ^^'^ ^ dXi dXj (3.11) 

J Pni \Xi jPm \Xj J 

. (3.12) 

Our method is to use In^^miXi^X^ as edge weights on a full graph and run Kruskal's 
algorithm until we encounter edges with negative weight. Let T be the set of all forests and 
'^n2{hj) — In2,m{Xi, Xj). The final forest is then 

= arg max {F) = arg min R^^ {pp) (3. 13) 

FeT FeT 

By building a forest on held-out data, we directly cross-vahdate over all forests. 
4. Statistical Properties 

In this section we present our theoretical results on risk consistency, structure selection 
consistency, and estimation consistency of the forest density estimate Pn = p-ik) ■ 

To establish some notation, we write a„ = f2(&n) if there exists a constant c such that 
On > cbn for sufficiently large n. We also write a„ x bn if there exists a constant c such that 
On < cbn and bn < can for sufficiently large n. Given a d-dimensional function / on the 
domain X, we denote its L2(P)-norm and sup-norm as 

\\f\\L2(P) = \ f P{x)dPx{x), ll/IU = sup \f{x)\ (4.1) 
y Jx xex 

where Px is the probability measure induced by X. Throughout this section, all constants 
are treated as generic values, and as a result they can change from line to line. 

In our use of a data splitting scheme, we always adopt equally sized splits for simplicity, 
so that rii — n2 — n/2, noting that this does not affect the final rate of convergence. 
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4.1' Assumptions on the density 

Fix /9 > 0. For any rf-tuple a = {ai, . . . , a^) G N'^ and x = (xi, . . . , Xa) G X, we define 
= n^=i ^ • Let denote the differential operator 



QOLi-\ had 



For any real- valued d-dimensional function / on A" that is [/3 J -times continuously differen- 



tiable at point Xq E X, let Pi'^] (x) be its Taylor polynomial of degree [/3J at point Xq: 



ai+-+ad<L/3J 

Fix L > 0, and denote by E(/3, L, r, xq) the set of functions / : A" — )■ R that are |_/3J-times 
continuously differentiable at xo and satisfy 

fi^)-Pf%i^) <L\\x-xo\\l, yxeB{xo,r) (4.4) 

where B{xo,r) = {x : \\x — xo\\2 < r} is the L2-ball of radius r centered at xq. The set 
E(/3, L, r, Xo) is called the L, r, a;o)-locally Holder class of functions. Given a set A, we 
define 

E(/3, L, r, A) = n,„eAE(^, L, r, Xq). (4.5) 
The following are the regularity assumptions we make on the true density function p*{x). 
Assumption 4.1. For any 1 < i < j < d, we assume 

(Dl) there exist Li > and L2 > such that for any c > the true bivariate and univariate 
densities satisfy 

p*{xi, Xj) e E (^f3, L2, c (log n/n)^ , Xi x Xj^ (4.6) 

and 

p*{xi) e E (/3,Li,c(logn/n)WT , A",) ; (4.7) 
(D2) there exists two constants ci and C2 such that 

Ci < inf p*{xi,Xj) < sup p*{xi,Xj) < C2 (4.8) 

Xi ,Xj EXiXXj 3.^. ^XiXXj 

/i- almost surely. 

These assumptions are mild, in the sense that instead of adding constraints on the joint 
density p*{x), we only add regularity conditions on the bivariate and univariate marginals. 
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4.. 2. Assumptions on the kernel 

An important ingredient in our analysis is an exponential concentration result for the kernel 
density estimate, due to Gine and Guillou (2002). We first specify the requirements on the 
kernel function K{-). 

Let (fi, A) be a measurable space and let J-" be a uniformly bounded collection of mea- 
surable functions. 

Definition 4.1. J-" is a bounded measurable VC class of functions with characteristics A 
and V if it is separable and for every probability measure P on (f2, A) and any < e < 1, 

iV(e||F|U,(P),^,l|-lU,(P)) < (7)', (4.9) 

where F(x) = supjgjr |/(x)| and A^(e,J-', || ■ ||l2(p)) denotes the e-covering number of the 
metric space {Q, \\ ■ \\l2{p))'i that is, the smallest number of balls of radius no larger than e 
(in the norm || ■ ||l2(p)) needed to cover J^. 

The one- dimensional density estimates are constructed using a kernel K, and the two- 
dimensional estimates are constructed using the product kernel 

K2{x,y) = K{x)-Kiy). (4.10) 

Assumption 4.2. The kernel K satisfies the following properties. 

//•oo 
K{u) du = 1, / K^{u) du < 00 and sup K{u) < c for some constant c. 
J-00 ueR 

(K2) K is a finite linear combination of functions g whose epigraphs epi{g) = {{s,u) : 

g[s) > u}, can be represented as a finite number of Boolean operations (union and 

intersection) among sets of the form {(s, u) : Q{s, u) > (f){u)}, where Q is a polynomial 

on M X M and is an arbitrary real function. 

(K3) K has a compact support and for any i > 1 and 1 < i' < [/3J 

j \tf \K{t)\dt < 00, and j \K{t)fdt < j t^' K{t)dt = 0. (4.11) 

Assumptions (Kl), (K2) and (K3) are mild. As pointed out by Nolan and Pollard (1987), 
both the pyramid (truncated or not) kernel and the boxcar kernel satisfy them. It follows 
from (K2) that the classes of functions 

J^i = j ( ) : M G M, /ii > }> (4.12 











\ hi 




(u-- 






\ h2 



J^2 = {i^K K {^-j—] ■.u,tER, h2>0} (4.13) 
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are bounded VC classes, in the sense of Definition 4.1. Assumption (K3) essentially says that 
the kernel K{-) should be P-valid; see Tsybakov (2008) and Definition 6.1 in Rigollet and Vert 
(2009) for further details about this assumption. 

We choose the bandwidths hi and h2 used in the one-dimensional and two-dimensional 
kernel density estimates to satisfy 

hi 
h2 

This choice of bandwidths ensures the optimal rate of convergence. 



logn 



1+2/3 



n 

\ogn\ 2W 

n 



(4.14) 
(4.15) 



4-3. Risk consistency 

Given the above assumptions, we first present a key lemma that establishes the rates of 
convergence of bivariate and univariate kernel density estimates in the sup norm. The proof 
of this and our other technical results are provided in Appendix A. 

Lemma 4.1. Under Assumptions 4-i one? 4-^! '^^^ choosing bandwidths satisfying (4.14) 
and (4.15), the bivariate and univariate kernel density estimates p{xi, Xj) andp{xk) in (3.3) 
and (3.4) satisfy 



max sup \p(xi, Xj) — p*(xi, Xi)\ = Op I \l 1 (4.16) 



and 



To describe the risk consistency result, let vjf" = Vd be the family of densities that are 
supported by forests with at most d — 1 edges, as already defined in (2.2). For Q < k < d — 1, 

Ik) 

we define V\ as the family of (i- dimensional densities that are supported by forests with at 
most k edges. Then 

VT ^Vf^ (4.18) 
Now, due to the nesting property (4.18), we have 

inf R{qF)> inf Riqp) > ■ ■ ■ > inf Rigp)- (4.19) 

We first analyze the forest density estimator obtained using a fixed number of edges k < d; 
specifically, consider stopping the Chow-Liu algorithm in Stage 1 after k iterations. This is 
in contrast to the algorithm described in 3.2, where the pruned tree size is automatically 
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determined on the held out data. While this is not very realistic in applications, since the 
tuning parameter k is generally hard to choose, the analysis in this case is simpler, and can 
be directly exploited to analyze the more complicated data-dependent method. 

Theorem 4.1 (Risk consistency). Letps,(k) be the forest density estimate with \E{f{''^)\ = k, 
obtained after the first k iterations of the Chow-Liu algorithm, for some k E {0, . . . , c? — 1}. 
Under Assumptions 4-i md 4-2, we have 

Note that this result allows the dimension d to increase at a rate o 

and the number of edges k to increase at a rate a ^a/ r2^/(i+^) / log , with the excess risk 
still decreasing to zero asymptotically. 

The above results can be used to prove a risk consistency result for the data-dependent 
pruning method using the data-splitting scheme described in Section 3.2. 

Theorem 4.2. Let p-(k) be the forest density estimate using the data- dependent pruning 

^d 

method in Section 3.2, and let pp(k) be the estimate with \E{Flf^)\ = k obtained after the 

d 

first k iterations of the Chow-Liu algorithm. Under Assumptions 4-1 and 4-2, we have 
where k* = argmino<fe<rf_i R(p^(k)). 

- - ^d 

The proof of this theorem is given in the appendix. A parallel result can be obtained for 
the method described in Section 3.3, which builds the forest by running Kruskal's algorithm 
on the heldout data. 

Theorem 4.3. Let Fn2 be the forest obtained using Kruskal's algorithm on held-out data, 
and let k = l-Fnal number of edges in F„2. Then 

where k* = \F*\ is the number of edges in the optimal forest F* = aigmmp^-p R{pf) . 
4-4- Structure selection consistency 

In this section, we provide conditions guaranteeing that the procedure is structure selection 
consistent. Again, we do not assume the true density p*{x) is consistent with a forest; rather, 
we are interested in comparing the estimated forest structure to the oracle forest which 
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minimizes the risk. In this way our result differs from that in Tan et al. (2009), although 
there are similarities in the analysis. 
By Proposition 2.2, we can define 

p*,k)=axgmmR{qF). (4.23) 

Thus f'^^ is the optimal forest within P^^'* that minimizes the negative log-likelihood loss. 
Let f'^'' be the estimated forest structure, fixing the number of edges at k; we want to study 
conditions under which 

P (Ff = ) ^ 1. (4.24) 

Let's first consider the population version of the algorithm — if the algorithm cannot re- 
cover the best forest Ff^ in this ideal case, there is no hope for stable recovery in the data 
version. The key observation is that the graph selected by the Chow-Liu algorithm only 
depends on the relative order of the edges with respect to mutual information, not on the 
specific mutual information values. Let 



S = |{(i,j),(fc,£)} : i < J and k < i, j and ij^kje {l,...,d}j. (4.25) 
The cardinality of £ is 

\£\ = 0{d^). (4.26) 

Let e = be an edge; the corresponding mutual information associated with e is denoted 
as Jg. If for all (e, e') G £, we have /g 7^ le', the population version of the Chow-Liu algorithm 

(k) 

will always obtain the unique solution F^ . However, this condition is, in a sense, both too 
weak and too strong. It is too weak because the sample estimates of the mutual information 
values will only approximate the population values, and could change the relative ordering 
of some edges. However, the assumption is too strong because, in fact, the relative order of 
many edge pairs might be changed without affecting the graph selected by the algorithm. 
For instance, when k > 2 and Jg and Jg' are the largest two mutual information values, it's 
guaranteed that e and e' will both be included in the learned forest F^'^^ whether /g > Jg/ or 
Jg < Jg/. 

Define the crucial set J <Z £ to be a set of pairs of edges (e, e') such that Jg 7^ Jg/ and 
flipping the relative order of Jg and Jg/ changes the learned forest structure in the population 
Chow-Liu algorithm, with positive probability. Here, we assume that the Chow- Liu algorithm 
randomly selects an edge when a tie occurs. 

The cardinality IJT"! of the crucial set is a function of the true density p*{x), and we can 
expect \ J\ <^ \£\. The next assumption provides a sufficient condition for the two-stage 
procedure to be structure selection consistent. 

Assumption 4.3. Let the crucial set J7 be defined as before. Suppose that 

min |J(X,;X,)-J(X,;X,)| >2L„ (4.27) 



16 



where L„ = Q 



'log 72 + logfi 



This assumption is satisfied in many cases. For example, in a graph with population 
mutual informations differing by a constant, the assumption holds. Assumption 4.3 is trivially 

satisfied if — )■ oo. However, if two pairs of edges belonging J' have the same 

log n + log a 

marginal distributions, the assumption may fail. 



Theorem 4.4 (Structure selection consistency). Let be the optimal forest within 
that minimizes the negative log-likelihood loss. Let f'^^ he the estimated forest with \Ep{k) \ = 

d 

k. Under Assumptions 4-1, 4-^> ^^i^d 4-3, we have 

P (^Ff ^ = ) j ^ 1 (4.28) 

as n oo. 

The proof shows that our method is structure selection consistent as long as the di- 
mension mcreases as d = o (exp(n^/(i+''))); in this case the error decreases at the rate 

o ^exp ^4 log d — c(log n) w log d j j . 



4-5. Estimation consistency 

Estimation consistency can be easily established using the structure selection consistency 
result above. Define the event Mk = {E^''^ = Fjf^^}. Theorem 4.4 shows that F{Ml) — )• as 
n goes to infinity. 

Lemma 4.2. Letp^(k) he the forest-based kernel density estimate for some fixed k G {0,...,(i— 
1}, and let 




(4.29) 



Under the assumptions of Theorem 4.4> 

= R{ppi.)) - Rip*,,,) (4.30) 

^d ^d ^d 

on the event Aik- 

Proof. According to Bach and Jordan (2003), for a given forest F and a target distribution 

p*(x), 

D{p*\\ qp) = D{p*\\ p*p) + D{p^ qp) (4.31) 
for all distributions qp that are supported by F. We further have 

D{p*\\q)= I p*{x)\ogp*{x) — I p* (x) log q{x)dx = / p* (x) log p* {x)dx + R{q) (4.32) 
Jx Jx Jx 
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for any distribution q. Using (4.31) and (4.32), and conditioning on the event M.k) "we have 
^(p;(.)l|P^w) = D{f\\p^,,,)-D{f\\p*,,,) (4.33) 
= I p*{x)\ogp*{x)dx + R{%(k)) - / p* (x) log p* {x)dx - R{p* ,k)) 

JX JX 

^d i'd 

which gives the desired result. □ 

The above lemma combined with Theorem 4.1 allows us to obtain the following estimation 
consistency result, the proof of which is omitted. 

Corollary 4.1 (Estimation consistency). Under Assumptions 4-1, 4-^> ^''^d 4-3, we have 

nf * II- ^ nil /logn + logc/ /logn + logrf\ 

^(P^(.)lbpw) = 0p \^k^ +d^ ^2/3/(1+2/.) j ■ (4-34) 

5. Tree Restricted Forests 

We now turn to the problem of estimating forests with restricted tree sizes. As discussed 
in the introduction, clustering problems motivate the goal of constructing forest structured 
density estimators where each connected component has a restricted number of edges. But 
estimating restricted tree size forests can also be useful in model selection for the purpose of 
risk minimization, since the maximum subtree size can be viewed as an additional complexity 
parameter. 

Definition 5.1. A t-restricted forest of a graph G is a subgraph Ft such that 

1. Ft is the disjoint union of connected components {Ti, ...,Tm}, each of which is a tree; 

2. |Tj| < t for each i < m, where |Tj| denotes the number of edges in the zth component. 

Given a weight We assigned to each edge of G, an optimal t-restricted forest Ft satisfies 

wiFt*) = locmx wiF) (5.1) 

where w{F) = 'Yle^p'^e is the weight of a forest F and J^tiG) denotes the collection of all 
t-restricted forests of G. 

For t = 1, the problem is maximum weighted matching. However, for t > 7, we show 
that finding an optimal t-restricted forest is an NP-hard problem; however, this problem 
appears not to have been previously studied. Our reduction is from Exact 3-Cover (X3C), 
shown to be NP-complete by Garey and Johnson (1979)). In X3C, we are given a set X, a 
family S of 3-element subsets of X, and we must choose a subfamily of disjoint 3-element 
subsets to cover X. Our reduction constructs a graph with special tree-shaped subgraphs 
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Algorithm 5.1 Approximate Max Weiglit t- Restricted Forest 
1: Input graph G with positive edge weights, and positive integer t > 2. 
2: Sort edges in decreasing order of weight. 
3: Greedily add edges in decreasing order of weight such that 

(a) the degree of any node is at most t + 1; 

(b) no cycles are formed. 

The resulting forest is F' = {Ti, T2, . . . , r„}. 
4: Output Ft = UjTreePartition(Tj , t). 



called gadgets, such that each gadget corresponds to a 3-element subset in S. We show that 
finding a maximum weight t-restricted forest on this graph would allow us to then recover a 
solution to X3C by analyzing how the optimal forest must partition each of the gadgets. 

Given the NP-hardness for finding optimal t- restricted forest, it is of interest to study 
approximation algorithms for the problem. Our first algorithm is Algorithm 5.1, which runs 
in two stages. In the first stage, a forest is greedily constructed in such a way that each 
node has degree no larger than t (a property that is satisfied by all t-restricted forests). 
However, the trees in the forest may have more than t edges; hence, in the second stage, 
each tree in the forest is partitioned in an optimal way by removing edges, resulting in a 
collection of trees, each of which has size at most t. The second stage employs a procedure 
we call TreePartition that takes a tree and returns the optimal t-restricted subforest. 
TreePartition is a divide-and-conquer procedure of Lukes (1974) that finds a carefully 
chosen set of forest partitions for each child subtree. It then merges these sets with the 
parent node one subtree at a time. The details of the TreePartition procedure are given 
in Appendix A. 

Theorem 5.1. Let Ft be the output of Algorithm 5.1, and let F* he the optimal t-restricted 
forest. Thenw{Ft) > ^w{F;). 

In Appendix A. 7, we present a proof of the above result. In that section, we also present 
an improved approximation algorithm, one based on solving linear programs, that finds a 
t-restricted forest F/ such that w{Fl) > ^w{F^). 

5.1. Pruning Based on t-Restricted Forests 

For a given t, after producing an approximate maximum weight t-restricted forest Ft using 
Vi, we prune away edges using V2. To do so, we first construct a new set of univariate and 
bivariate kernel density estimates using V2, as before, Pn2{xi) and Pn2{xi,Xj). Recall that we 
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Algorithm 5.2 t- Restricted Forest Density Estimation 



1: Divide data into two halves Vi and X'2- 

2: Compute kernel density estimators and Pn2 for all pairs and single variable marginals. 

3: For all pairs compute In^iXi, Xj) according to (3.5) and In2,ni{Xi, Xj) according to (5.3). 

4: For t = 0, . . . , tfinai where ifinai is chosen based on the application 

1. Compute or approximate (for t >2) the optimal i-restricted forest Ft using as edge weights. 

2. Prune Ft to eliminate all edges with negative weights In2,ni- 

5: Among all pruned forests ppt, select t = argminQ^f^j^.^^^^ Rn2iPpJ- 



define the "cross-entropies" of the kernel density estimates for each pair of variables as 

Tn2,ni{Xi,Xj) = I Pn^ {Xi , Xj ) log ^ ^^1' ^ dXi dXj (5.2) 
J Pni{Xi)Pni{Xj ) 

^ , ) log ^ -^"^ . (5.3) 

We then eliminate all edges in Ft for which /„2 „j(Xj, Xj) < 0. For notational simplicity, 
we denote the resulting pruned forest again by Ft. 

To estimate the risk, we simply use Rn2(j>pJ as defined in (3.6), and select the forest F-j- 
according to 

t = arg min R^^ (pp ). (5.4) 

0<t<d~l 

The resulting procedure is summarized in Algorithm 5.2. 

Using the approximation guarantee and our previous analysis, we have that the population 
weights of the approximate t-restricted forest and the optimal forest satisfy the following 
inequality. We state the result for a general c-approximation algorithm; for the algorithm 
given above, c = 4, but tighter approximations are possible. 

Theorem 5.2. Assume the conditions of Theorem 4-1- Fort > 2, let Ft be the forest con- 
structed using a c-approximation algorithm, and let F* be the optimal forest; both constructed 
with respect to finite sample edge weights Wm = Im ■ Then 

MP,)>i^F;) + o,((r + S)^!^i^) (5.5) 

where k and k* are the number of edges in Ft and F^ , respectively, and w denotes the 
population weights, given by the mutual information. 

As seen below, although the approximation algorithm has weaker theoretical guarantees, 
it out-performs other approaches in experiments. 
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6. Experimental Results 

In this section, we report numerical results on both synthetic datasets and microarray data. 
We mainly compare the forest density estimator with sparse Gaussian graphical models, 
fitting a multivariate Gaussian with a sparse inverse covariance matrix. The sparse Gaussian 
models are estimated using the graphical lasso algorithm (glasso) of Friedman, Hastie and Tibshirani 

(2007) , which is a refined version of an algorithm first derived by Banerjee, El Ghaoui and d'Aspremont 

(2008) . Since the glasso typically results in a large parameter bias as a consequence of the 
ii regularization, we also compare with a method that we call the refit glasso, which is a 
two-step procedure — in the first step, a sparse inverse covariance matrix is obtained by the 
glasso; in the second step, a Gaussian model is refit without ii regularization, but enforcing 
the sparsity pattern obtained in the first step. 

To quantitatively compare the performance of these estimators, we calculate the log- 
likelihood of all methods on a held-out dataset V2. With and flm denoting the estimates 
from the Gaussian model, the held-out log-likehhood can be explicitly evaluated as 

W = -^^11 [l^^^'^ - l^nyn^AX^^^ - + \ log (^^^ I . (6.1) 

For a given tree structure F, the held-out log-likelihood for the forest density estimator is 

where PnA') the corresponding kernel density estimates using the plug-in bandwidths. 

Since the held-out log-likelihood of the forest density estimator is indexed by the number 
of edges included in the tree, while the held-out log-likelihoods of the glasso and the refit 
glasso are indexed by a continuously varying regularization parameter, we need to find a 
way to calibrate them. To address this issue, we plot the held-out log-likelihood of the forest 
density estimator as a step function indexed by the tree size. We then run the full path of the 
glasso and discretize it according to the corresponding sparsity level, i.e., how many edges 
are selected for each value of the regularization parameter. The size of the forest density 
estimator and the sparsity level of the glasso (and the refit glasso) can then be aligned for a 
fair comparison. 

6.1. Synthetic data 

We use a procedure to generate high dimensional Gaussian and non-Gaussian data which 
are consistent with an undirected graph. We generate high dimensional graphs that contain 
cycles, and so are not forests. In dimension d = 100, we sample ni = ^2 = 400 data 
points from a multivariate Gaussian distribution with mean vector fi = (0.5, . . . , 0.5) and 
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inverse covariance matrix f2. The diagonal elements of f2 are all 62. We then randomly 
generate many connected subgraphs containing no more than eight nodes each, and set the 
corresponding non-diagonal elements in Q at random, drawing values uniformly from —30 
to —10. To obtain non-Gaussian data, we simply transform each dimension of the data by 
its empirical distribution function; such a transformation preserves the graph structure but 
the joint distribution is no longer Gaussian (see Liu, Lafferty and Wasserman (2009)). 

To calculate the pairwise mutual information I{Xi; Xj), we need to numerically evaluate 
two-dimensional integrals. We first rescale the data into [0, l]'' and calculate the kernel density 
estimates on a grid of points; we choose m = 128 evaluation points x^p < xf^ < ■ ■ ■ < x'f^^ for 
each dimension z, and then evaluate the bivariate and the univariate kernel density estimates 
on this grid. 

There are three different kernel density estimates that we use — the bivariate kde, the 
univariate kde, and the marginalized bivariate kde. Specifically, the bivariate kernel density 
estimate on Xi.Xj based on the observations {x'f\x'^^''}s&Vi is defined as 

using a product kernel. The bandwidths h2i, are chosen as 

, 1 nc • Q'fc,0.75 — Q'A:,0.25 1 -1/(2/3+2) /rA\ 

/i2fc = 1-06 ■ mm <^ (Tfc, (''^ /^P+)^ (6.4) 

where is the sample standard deviation of {xjf^}s^Vi and ^fc,o.25 are the 75% and 

25% sample quantiles of {X^''''}<ier>i- 

In all the experiments, we set /3 = 2, such a choice of f3 and the "plug-in" bandwidth h2k 
(and hik in the following) is a very common practice in nonparametric Statistics. For more 
details, see Fan and Gijbels (1996) and Tsybakov (2008). 

Given an evaluation point Xk, the univariate kernel density estimate p{xk) based on the 
observations {xj:^}s(zDi is defined as 




where hik > is defined as 

, 1 -1^ ?fc,0.75 ~ 9fc,0.25 1 _i/(2fl+l) /r r\ 

/iifc = 1.06 ■ mm <^ cTfc, ^ — >-n p+ K (6.6) 

Finally, the marginal univariate kernel density estimate PM{xk) based on the observations 

{xjf^}s(zVi is defined by integrating the irrelevant dimension out of the bivariate kernel 
density estimates p{xj,Xk) on the unit square [0, 1]^. Thus, 

PM{xk) = y^Mxf ,Xk). (6.7) 
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With the above definitions of the bivariate and univariate kernel density estimates, we con- 
sider estimating the mutual information /(Xj; Xj) in three different ways, depending on 
which estimates for the univariate densities are employed. 



^ mm 

k'=ie'=i 

— - ^^(xf' ^) logp(xf' ^) 7 y]p(a;f ^) logp(5 

— 1 ^-^ m — 1 ^-^ ■' 



fc'=i i'=i 



/medium(Xi,X,) = — — ^ ^ p(xf \ xf ^ ) log — ^^;y-^-^ . (6.9) 

^ mm 

/siow(X„X,) = 7^^^3^EE^(^r\4''^)logfe^''^<^) - (6-10) 

k'=ie'=i 

TTl TTl 

m — 1 ^-^ m — I ^-^ ■' ■' 

k'=l £'=1 

The terms "fast," "medium" and "slow" refer to the theoretical statistical rates of conver- 
gence of the estimators. The "fast" estimate uses one- dimensional univariate kernel density 
estimators wherever possible. The "medium" estimate uses the one-dimensional kernel den- 
sity estimates in the denominator of p{xi,Xj)/{p{xi)p{xj), but averages with respect to the 
bivariate density. Finally, the "slow" estimate marginalizes the bivariate densities to esti- 
mate the univariate densities. While the rate of convergence is the two-dimensional rate, the 
"slow" estimate ensures the consistency of the bivariate and univariate densities. 

Figure 1 compares /fast, ^medium, and Jgiow on different pairs of variables. The boxplots 
are based on 100 trials. Compared to the ground truth, which can be computed exactly in 
the Gaussian case, we see that the performance of /medium and /glow is better than that of 
/fast- This is due to the fact that simply replacing the population density with a "plug-in" 
version can lead to biased estimates; in fact, /fast is not even guaranteed to be non- negative. 
In what follows, we employ /medium for all the calculations, due to its ease of computation and 
good finite sample performance. Figure 2 compares the bivariate fits of the kernel density 
estimates and the Gaussian models over four edges. For the Gaussian fits of each edge, we 
directly calculate the bivariate sample covariance and sample mean and plug them into the 
bivariate Gaussian density function. From the perspective and contour plots, we see that the 
bivariate kernel density estimates provide reasonable fits for these bivariate components. 

A typical run showing the held-out log-likelihood and estimated graphs is provided in Fig- 
ure 3. We see that for the Gaussian data, the refit glasso has a higher held-out log-likelihood 
than the forest density estimator and the glasso. This is expected, since the Gaussian model 
is correct. For very sparse models, however, the performance of the glasso is worse than that 
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Fig 1: (Gaussian example) Boxplots of /fast; -^medium) and /slow on three different pairs of variables. 
The red-dashed horizontal lines represent the population values. 
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Fig 2: Perspective and contour plots of the bivariate Gaussian fits vs. the kernel density 
estimates for two edges of a Gaussian graphical model. 
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Fig 3: Synthetic data. Top-left Gaussian, and top-right non-Gaussian: Held-out log-hkehhood 
plots of the forest density estimator (black step function), glasso (red stars), and refit glasso 
(blue circles), the vertical dashed red line indicates the size of the true graph. Bottom plots 
show the true and estimated graphs for the Gaussian (second row) and non-Gaussian data 
(third row). 



25 



of the forest density estimator, due to the large parameter bias resulting from the ii regu- 
larization. We also observe an efficiency loss in the nonparametric forest density estimator, 
compared to the refit glasso. The graphs are automatically selected using the held-out log- 
likelihood, and we see that the nonparametric forest-based kernel density estimator tends 
to select a sparser model, while the parametric Gaussian models tend to overselect. This 
observation is new and is quite typical in our simulations. Another observation is that the 
held-out log-likelihood curve of the glasso becomes fiat for less sparse models but never goes 
down. This suggests that the held-out log-likelihood is not a good model selection criterion 
for the glasso. For the non-Gaussian data, even though the refit glasso results in a reasonable 
graph, the forest density estimator performs much better in terms of held-out log-likelihood 
risk and graph estimation accuracy. 

To compare with t-restricted forests, we generated additional Gaussian and non-Gaussian 
synthetic data as before except on a different graph structure. In Figure 4, we use 400 
training examples while varying the size of heldout data to compare the log-likelihoods of 
four different methods; the log-likelihood is evaluated on a third large dataset. In Figure 5, we 
consider only non-Gaussian data, use 400 training data and 400 heldout data, and generate 
graphs with best heldout log-likelihood across the four methods. We compute bandwidth, 
heldout log-likelihood, and mutual information same as before. 
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Fig 4: Log-likelihood comparison of various methods: (left white) MST on Training Data 
with Pruning (gray) MST on Heldout Data (black) t-Restricted Graph (blue) Refit Glasso 



We observe that although creating a maximum spanning tree (MST) on the held-out data 
is asymptotically optimal; it can perform quite poorly. Unless there are copious amount of 
heldout data, held-out MST overfits on the heldout data and tend to give large graphs; in 
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(d) (e) 

Fig 5: Graphs generated on non-Gaussian Data: (a) True Graph, (b) t-Restricted Forest (c) 
MST on Heldout Data (d) MST on Training Data with Pruning (e) Refit Glasso 

contrast, t- restricted forest has the weakest theoretical guarantee but it gives the best log- 
hkelihood and produces sparser graphs. It is not surprising to note that MST on heldout 
data improves as heldout data size increases. Somewhat surprisingly though, Training-MST- 
with-pruning and t-restricted forest appear to be insensitive to the heldout data size. 

6.2. Microarray data 

6.2.1. Arabidopsis thaliana Data 

In this example, we consider a dataset based on Affymetrix GeneChip microarrays for the 
plant Arabidopsis thaliana, (Wille et ai, 2004). The sample size is n = 118. The expression 
levels for each chip are pre-processed by a log-transformation and standardization. A subset 
of 40 genes from the isoprenoid pathway are chosen, and we study the associations among 
them using the glasso, the refit glasso, and the tree-based kernel density estimator. 
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Fig 6: Results on microarray data. Top: held-out log-likelihood (left) and its zoom-in (right) 
of the tree-based kernel density estimator (black step function), glasso (red stars), and refit 
glasso (blue circles). Bottom: estimated graphs using the tree-based estimator (left) and 
glasso (right). 
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From the held-out log-likelihood curves in Figure 6, we see that the tree-based kernel 
density estimator has a better generalization performance than the glasso and the refit glasso. 
This is not surprising, given that the true distribution of the data is not Gaussian. Another 
observation is that for the tree-based kernel density estimator, the held-out log-likelihood 
curve achieves a maximum when there are only 35 edges in the model. In contrast, the 
held-out log-likelihood curves of the glasso and refit glasso achieve maxima when there are 
around 280 edges and 100 edges respectively, while their predictive estimates are still inferior 
to those of the tree-based kernel density estimator. 

Figure 6 also shows the estimated graphs for the tree-based kernel density estimator and 
the glasso. The graphs are automatically selected based on held-out log-likehhood. The two 
graphs are clearly different; it appears that the nonparametric tree-based kernel density esti- 
mator has the potential to provide different biological insights than the parametric Gaussian 
graphical model. 

6.2.2. HapMap Data 

This dataset comes from Nayak et al. (2009). The dataset contains Affymetrics chip mea- 
sured expression levels of 4238 genes for 295 normal subjects in the Centre d'Etude du Poly- 
morphisme Humain (CEPH) and the International HapMap collections. The 295 subjects 
come from four different groups: 148 unrelated grandparents in the CEPH-Utah pedigrees, 
43 Han Chinese in Beijing, 44 Japanese in Tokyo, and 60 Yoruba in Ibadan, Nigeria. Since 
we want to find common network patterns across different groups of subjects, we pooled the 
data together into a n = 295 by = 4238 numerical matrix. 

We estimate the full 4238 node graph using both the forest density estimator (described 
in Section 3.1 and 3.2) and the Meinshausen-Biihlmann neighborhood search method as 
proposed in (Meinshausen and Biihlmann, 2006) with regularization parameter chosen to 
give it about same number as edges as the forest graph. 

To construct the kernel density estimates p{xi,Xj) we use an array of Nvidia graphical 
processing units (GPU) to parallelize the computation over the pairs of variables Xj and 
Xj. We discretize the domain of {Xi,Xj) into a 128 x 128 grid, and correspondingly employ 
128 X 128 parallel cells in the GPU array, taking advantage of shared memory in CUDA. 
Parallelizing in this way increases the total performance by approximately a factor of 40, 
allowing the experiment to complete in a day. 

The forest density estimated graph reveals one strongly connected component of more 
than 3000 genes and various isolated genes; this is consistent with the analysis in Nayak et al. 
(2009) and is realistic for the regulatory system of humans. The Gaussian graph contains 
similar component structure, but the set of edges differs significantly. We also ran the t- 
restricted forest algorithm for t = 2000 and it successfully separates the giant component 
into three smaller components. For visualization purposes, in Figure 7, we show only a 934 
gene subgraph of the strongly connected component among the full 4238 node graphs we 
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estimated. More detailed analysis of the biological implications of this work will left as a 
future study. 




Fig 7: A 934 gene subgraph of the full estimated 4238 gene network. Left: estimated forest 
graph. Right: estimated Gaussian graph. Red edges in the forest graph are missing from the 
Gaussian graph and vice versa; the blue edges are shared by both graphs. Note that the 
layout of the genes is the same for both graphs. 



7. Conclusion 

We have studied forest density estimation for high dimensional data. Forest density estima- 
tion skirts the curse of dimensionality by restricting to undirected graphs without cycles, 
while allowing fully nonparametric marginal densities. The method is computationally sim- 
ple, and the optimal size of the forest can be robustly selected by a data-splitting scheme. 
We have established oracle properties and rates of convergence for function estimation in 
this setting. Our experimental results compared the forest density estimator to the sparse 
Gaussian graphical model in terms of both predictive risk and the qualitative properties of 
the estimated graphs for human gene expression array data. Together, these results indicate 
that forest density estimation can be a useful tool for relaxing the normality assumption in 
graphical modeling. 
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Appendix A: Proofs 
A.l. Proof of Lemma 4 • 1 

We only need to consider tlie more complicated bivariate case (4.16); the result in (4.17) 
follows from the same line of proof. First, given the assumptions, the following lemma can 
be obtained by an application of Corollary 2.2 of Gine and Guillou (2002). For a detailed 
proof, see Rinaldo and Wasserman (2009). 

Lemma A.l. (Gine and Guillou, 2002) Let p he a bivariate kernel density estimate using 
a kernel K{-) for which Assumption 4-2 holds and suppose that 

sup sup / K2{u)p*{t — uh2)du < D < oo. (A.l) 

1. Let the bandwidth h2 be fixed. Then there exit constants L > and C > 0, which 
depend only on the VC characteristics of in (4.13), such that for any ci > C and 
< e < CiD /\\K2\\oo! there exists no > which depends on e, D, ||-ft'2||oo ond the VG 
characteristics of K2, such that for all n> Uq, 

P f sup \piu) - Epiu)\ > 2e] < Lexp ( _ 1 Ml + Ml^O) !^ 1 ^^.2) 



L ci D 



2. Let /i2 — ^ m such a way that n/i2/log/i2 00, and let e ^ so that 



where Vn = Vl{h2^). Then (A. 2) holds for sufficiently large n. 

From (-D2) in Assumption 4.1 and {Kl) in Assumption 4.2, it's easy to see that (A.l) is 
satisfied. Also, since 

(A.4) 

n J 

it's clear that nh\/\ogh2 — )• 00. Part 2 of Lemma A.l shows that there exist C2 and C3 such 
that 



P 



sup \p{xi,Xj) —Wf){xi.,Xj)\ > - ) < C2exp (—czn^+f{\ogn)^+i^e^\ (A.5) 



for all e satisfying (A. 3). 

This shows that for any i,j G {1, . . . ,d} with i j, the bivariate kernel density estimate 
p{xi,Xj) is uniformly close to Ep(xi,Xj). Note that Kp{xi,Xj) can be written as 

Xi, Xj) = j -^K ^ K (^J^-J^^ P*iui, Vj) duidvj. (A. 6) 
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The next lemma, from RigoUet and Vert (2009), provides a uniform deviation bound on 
the bias term Ep(xj, Xj) — p*{xi, Xj). 

Lemma A. 2. (Rigollet and Vert, 2009) Under (-D1) in Assumption 4-1 and {K3) in As- 
sumption 4-2, we have 



sup \Ep{xi,Xj) -p*{x^,Xj)\< Lih^ {u^ + vY^^K 



where L is defined in (-D1) of Assumption 4-1- 



{u)K{v) dudv. (A.7) 



Let C4 = Li {u + V K{u)K{v) dudv. From the discussion of Example 6.1 in 

Rigollet and Vert (2009) and {Kl) in Assumption 4.2, we know that C4 < 00 and only 
depends on K and /3. Therefore 

P ( sup \p\xi, Xj) - ¥${xi, x,)| > ^ ) = (A.8) 

{xi,Xj)&Xi xXj 

for e > 4c4/i2- 

The desired result in Lemma 4.1 is an exponential probability inequality showing that 
p{xi,Xj) is close to p*{xi,Xj). To obtain this, we use a union bound: 

P max sup \pixi,Xj) ~ p*{xi,Xj)\ > e 

y(j,j)G{l,...,d}x{l,...,d} (xi,Xj)eXiXXj 

< d'^F [ sup \p{xi,Xj) — Ep(xj,Xj)| > - 

\^{xi,Xj)£XiXXj 

+ d^F I sup \p*{xi,Xj) - Ep{xi,Xj)\ > I • (A. 9) 

\{xi,Xj)GXiXXj 

Choosing 



the result directly follows by combining (A. 5) and (A.8) 
A. 2. Proof of Theorem 4-i 

First, from (-D2) in Assumption 4.1 and Lemma 4.1, we have for any i 7^ j, 



p{xi,Xj) ^ / /logn + logrf 



max sup ^ , " ■'\ -1] =0p\ \ ^ . (A.ll) 

{^,j)e{l,...,d}^{l,...,d}(x,,x,)eX,xX,\P*ix^,Xj) J \V n^^A/^+l) ' 
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The next lemma bounds the deviation of R{pf) from R{p*p) over different choices of 
F e with \E{F)\ < k. In the following, we let 



J-f = {FG J-,:|i?(F)|<A;} 

denote the family of d-node forests with no more than k edges. 
Lemma A. 3. Under the assumptions of Theorem 4-1, we have 



(A. 12) 



\ Ml n \ u i^ogn + logd \ogn + \ogd 

sup ^ \R{pj.) - R{pp)\ =Op\^k^ + ^2/3/(1+2/.) 



(A 13) 



Proof. For any F G we have 



\Ripp) - R{p*p) 
< 



{i,j)£E{F) \^^iX^i 



p(^X{^ Xj 



Ai 



f / p*{xk)\ogp*{xk)dxk- I p{xk) log p{xk)dxk 



A2 



Defining p*^ = p*{xi, Xj) and pij = p{xi, xj), we further have 
A, = 

\{i,j)€E{F) 



(A.14) 



(A15) 



sup \p*{xi,Xj) — p{xi,Xj)\ / log p{xi,Xj)dxidxj 



where we use the fact the univariate and bivariate densities are assumed to have the same 
smoothness parameter (3, therefore the univariate terms are of higher order, and so can be 
safely ignored. 

It now suffices to show that 



Ai=Op\k 



log n + log d 



(A.16) 



and 



Ao = Op\d 



log n + log d 

^2/3/(1+2/3) 



(A.17) 
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In the sequel, we only prove (A. 16); (A. 17) follows in the same way. 

To show (A. 16), using the fact that niax(jj)g{i _rf}x{i,...,rf} Xj) < C2, it's sufficient to 
prove that 



and 



loff n + log d 



(i,j)G{i,...,d|x{i,...,d| - \v n/'/('3+i) '■ (A.19) 

Equation (A. 18) directly follows from (4.16) in Lemma 4.1, while (A.19) follows from the 
fact that, for any densities p and q, where q is strictly positive, 

D{p\\q)= [ ^log^q{x)dx. 
J q{x) q{x) 

By a Taylor expansion, for x ~ 1, 

X logx = (x — 1) + o (x — 1) 

and we then have 



(A.20) 
(A.21) 



DiPij \\Pij) = Op I sup \p*ixi, Xj) - p{xi, Xj) I I . 

{xi,Xj)£XiXXj 



(A.22) 
□ 



The desired result follows by combining (A. 18) and (A.19). 

The next auxiliary lemma is also needed to obtain the main result. It shows that R{pf) 
does not deviate much from R{pf) uniformly over different choices of F e J^jf'^ ■ 

Lemma A. 4. Under the assumptions of Theorem 4-1, we have 



wf- \ dY- M n 17 i^ogn + hgd /logn + logci 

sup ^ \RipF) - R{pf) I =Op^k^ + ^2/3/(1+2/.) 



(A.23) 



Proof. Following the same line of argument as in Lemma A. 3, we have for all F e 
\R{Pf)-R{Pf)\ (A.24) 



< 



+ 



diOC'idiOC'i 



f / p*(xfc) logp(xfc)(ixfc - / p(xfc) logp(xfc)(ixfc 



Op I E 

M,j)eE{F) 



sup \p* {xi ,Xj) — p{xi ,Xj)\ I log p{xi , Xj ) dxidxj 



^[ P* (xk) log p{xk)dxk- p{xk) log p{xk)dxk 

k&Vp ^ '^'^ 
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From (A. 11), we get that 

max log < max{| logC2|, I logcil} + 1 (A. 25) 

(i,i)6{l,...,d}x{l,...,d} 

for large enough n. The result then directly follows from (4.16) and (4.17) in Lemma 4.1. □ 

The proof of the main theorem follows by repeatedly applying the previous two lemmas. 
As in Proposition 2.2, with 

p* (fe) = argmini?(gF), (A.26) 



we have 



R{p^w ) - RiPpW ) + RiPpw ) - R{pl(k) ) (A.27) 



d d d 



T}{- \ T}( * \ ^ n \ 1 ^ogn + \ogd /logra + logrf , 
i?(Ppw)-^(p,w) + Op ( k\l ^^^^^^^^ ^2/3/(1+2/3) 1 (A-28) 



^ T?^-- \ u( * \ ^ n { 1 ^ogn + hgd \ogn + \ogd . 

< Rip^i.))-Rip^,,,) + Op\k\l ^^^^^^^^ +dd I (A.29) 



* \ * \ I n I u ^ogn + \ogd \ogn + \ogd , 

- +Op\^k^ ^/3/(/3+l) + ^2/3/(1+2/3) 1 (A-30) 



I , /logn + log(i /logn + logd , 

I ^ V ^/3/(/3+l) + ^2/3/(1+2/3) 1 • (A-31) 



- d 



where (A. 28) follows from Lemma A. 4, (A.29) follows from the fact that p^{k) is the minimizer 
of -R(-), and (A. 30) follows from Lemma A. 3. 

A. 3. Proof of Theorem 4.2 

To simplify notation, we denote 



<Pn{k) = kJ^-^j—^ (A.32) 



/loi 


2; n + log 




n/3/(/3+i) 


'loj 


5 + log 


■.d 



^«(^) = ^2/3/(1+2/3) • (A-33) 

Following the same proof as Lemma A. 4, we obtain the following. 
Lemma A. 5. Under the assumptions of Theorem ^.1, we have 

sup \R{pF)-Rn,{PF)\=Op((Pn{k)+^M^- (A.34) 
where is the held out risk. 
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To prove Theorem 4.2, we now have 

R{PpCk)) - RiPpi^*)) = R{Pp(k)) - Rn^iPpCk)) + Rn^iPpw) - R{Pp(kn) (A.35) 



- d d d 



= Op{Mk)+Md)) + RnAPBw)-RiPF(^*)) (A.36) 

^d 

< Op{(f)n(k)+^n{d)) + Rn,{PB(^'))-R{PB('^n) (A.37) 

^d 

= Op(^Mk) + Mk*)+Md)). (A.38) 
where (A.37) follows from the fact that k is the minimizer of i?„2(-). 

A. 4- Proof of Theorem 4-3 

Using the shorthand 



Mk) = kd t./n^m (A.39) 



'lo| 


g n + log 


:d 




'lO| 


5 n + log 


:d 



^nid) = d^ (A.40) 
We have that 

R{ppJ-Ripp.) = R{ppJ - Rn.ippJ + Rn.ippJ - R{PF') (A.41) 
= OpiMk)+Md)) + Rn,ippJ-RiPF*) (A.42) 
< Op{Mk)+Md)) + Rn,{pF*)-R{pF*) (A.43) 
= Op{<f)n(k) + Mkl+Md)) (A.44) 

(A.45) 

where line A.43 follows because is the minimizer of i?„2(')- 
A. 5. Proof of Theorem 4-4 

We begin by showing an exponential probability inequality on the difference between the 
empirical and population mutual informations. 

Lemma A. 6. Under Assumptions 4-1, 4-^> there exist generic constants C5 and cq satisfying 
P {\I{X,;Xj) - T{X,;X,)\ > e) < C5 exp (-Cgn^ (log n)We2j . (A.46) 
for arbitrary i,j G {1, . . . , d} with i 7^ j , and e so that 



where rn = VL{h2 
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Proof. For any e = i7 I -t / — ^pp- | , we have 



v(\I(X,;X,)-T{Xf,X,)\>e^ 

+ P(| / {f{xi,x^)\ogf{x^f{xj)-^{xi,Xj)\o<gy{x^^{x^)^ (A.4J 

Since the second term of (A. 48) only involves univariate kernel density estimates, this term 
is dominated by the first term, and we only need to analyze 



P ( I / (p*(xi,Xj)logp*(xj,Xj) -p(xi,Xj)logp(xi,Xj))(ixi(iXj| > ^ ) • (A.49) 



The desired result then follows from the same analysis as in Lemma A. 3. □ 
Let 

^" - " [^^) (A.50) 

be defined as in Assumption 4.3. To prove the main theorem, we see the event ^ F^^^ 
implies that there must be at least exist two pairs of edges and {k,£), such that 

sign(/(X„ X,) - I{Xk, X,)) ^ sign (/(X„ X,) - T{Xk, X,)) . (A.51) 

Therefore, we have 



p (^F^ + FT) (A.52) 
< P ((/(X„X,) -/(X,,X,)) ■ (/(X„X,) -/(X,,X,)) < 0, for some (A:,^)) • 

With d nodes, there can be no more than pairs of edges; thus, applying a union bound 
yields 



P ((/(X„X,) - /(Xfe,X,)) ■ (/(X„X,) - /(X,,X,)) < 0, for some (z, j), (fc,^)) (A.53) 

— max 

2 m)AKi))^j 



< - max P ((/(X„X,) - /(Xfc,X,)) • (/(X„X,) - /(Xfc,X,)) < o) . (A.54) 
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Assumption 4.3 specifies that 

mill - IiXk,Xe)\ > 2L„. (A.55) 

i{i,j),ik,e))ej 

Therefore, in order for (A. 51) hold, there must exist an edge G JT" such that 

\I{Xi,Xj)-T{X„Xj)\>L^. (A.56) 

Thus, we have 

P ((/(X„X,) - I{X,,Xe)) ■ (/(X„X,) - T{Xk,Xe)) < o) (A.57) 

< max P (\I{Xi, Xj) - /(Xi, X,)| > L„) (A.58) 

< c^exp (^-CQn^(\ogn)^Llj . (A. 59) 



max 



where (A. 59) follows from Lemma A. 6. 

Chaining together the above arguments, we obtain 



p (^pW _^ j (A.60) 

< P ((/(X„X,) - /(Xfc,X,)) ■ (T{X,,Xj) - T{Xk,Xe)) < 0, for some (fc,£)) 

< - max P ((l{X„X,) - I{Xk,Xe)) ■ (T{X,,X,) - /(Xfc,X,)) < o) (A.61) 

< rf^ max ¥ (\I{X„Xj) -T{Xi,Xj)\> Ln) (A.62) 

< (i^csexp (^-C6nW(logn)WL^j (A. 63) 

= o ^csexp ^41og(i - C6(logn)W logdj j (A. 64) 

= o(l). (A.65) 

The conclusion of the theorem now directly follows. 

A. 6. Proof of NP-hardness oft-Restricted Forest 

We will reduce an instance of exact 3-cover (X3C) to an instance of finding a maximum 
weight t- restricted forest (t-RF). 

Recall that in X3C, we are given a finite set X with |X| = 3q and a family of 3-element 
subsets of X, iS = {S* C X : \S\ = 3}. The objective is to find a subfamily S' C S such 
that every element of X occurs in exactly one member of S' , or to determine that no such 
subfamily exists. 

Suppose then we are given X = {xi, . . . , Xn} and 5 = {S* C X : l^l = 3}, with m = \S\. 
We construct the graph G in an instance of t-KF as follows, and as illustrated in Figure 8. 
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other gadgets 





Fig 8: Gadget constructed in the reduction from X3C 



For each x G X, add an element node to G. For each S* G 5, construct a gadget, which 
is a subgraph comprised of a nexus node, three junction nodes, and three lure nodes; see 
Figure 8. We assign weights to the edges in a gadget in the following manner: 

w (element, junction) = 2 

w (nexus, lurei) = 5 

iy(lurei, lure2) = 10 

w(lure2, lures) = 10 

w (nexus, junction) = > 31m. 

Note that the weight is chosen to be strictly greater than the weight all of the non-nexus- 
junction edges in the graph combined. To complete the instance of t-RF, let t = 7. 

Lemma A. 7. Suppose G is a graph constructed in the transformation from X3C described 
above. Then must contain all the nexus-junction edges. 

Proof. The set of all nexus-j unction edges together form a well-defined t-restricted forest, 
since each subtree has a nexus node and 3 junction nodes. Call this forest F. If some forest 
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F' is missing a nexus-j unction edge, then F' must have weight strictly less than F, since 
is larger than the sum of all of the non- nexus-j unction edges. □ 

Lemma A. 8. Each subtree in F^ can contain at most one nexus node. 

Proof. Suppose a subtree T in F* contains two nexus nodes. Then it must contain 6 junc- 
tion nodes by Lemma A. 7. Thus, T contains at least 8 nodes, and therefore violates the 
t- restriction constraint. □ 

Lemma A. 9. For each nexus node contained in F^ , the corresponding three junction nodes 
are either connected to all or none of the three neighboring element nodes. 

Proof. By the previous two Lemmas A. 7 and A. 8, each subtree is associated with at most one 
gadget, and hence at most one S" G 5, and moreover each gadget has as least one associated 
subtree. 

Without loss of generality, we consider a region of the graph local to some arbitrary 
subtree. By the size constraint, a subtree cannot contain all the adjacent element nodes and 
all the lure nodes. 

We now perform analysis: 

1. If a subtree contains no element nodes and all the lure nodes, then it has weight 3A+25. 
Call this an off configuration. 

2. If a subtree contains two element nodes, and a second subtree of three nodes contains 
all the lure nodes, then the total weight of both subtrees is 3 A -|- 24. This is suboptimal 
because we can convert to an off configuration and gain additional weight without 
affecting any other subtrees. Hence, such a configuration cannot exist in F^. 

3. If a subtree contains two element nodes and lurei, and a second subtree contains just 
lure2 and lures, then the total weight of the two subtrees is 3A + 19. This is again 
suboptimal. 

4. If a subtree contains an element node and both lurei and lure2, then there cannot 
be a second subtree in region local to the gadget. The weight of this one subtree is 
(3 A + 2 + 5 + 1 0) = 3A + 17, which is suboptimal. 

5. If a subtree contains all three element nodes and no lure nodes, and a second subtree 
contains all the lure nodes, then the total weight is (3A + 6) + 20 = 3A + 26. Call this 
an ON configuration. 

Thus, we see that each gadget in F^ must be either an ON or an off configuration. □ 

Recall that each gadget corresponds to a 3-element subset S in the family S. Since a 
gadget in an ON configuration has greater weight than a gadget in an off configuration, an 
optimal t-RF will have as many gadgets in the ON configuration as possible. Thus, to solve 
X3C we can find the optimal t-RF and, to obtain a subcover 5', we place all S into S' that 
correspond to ON gadgets in the forest. By Lemma A. 8 each subtree can contain at most one 
nexus node, which implies that each ON gadget is connected to element nodes that are not 
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connected to any other ON gadgets. Thus, this results in a subcover for which each element 
of X appears in at most one S & S'. 

A . 7. Proof of Theorem 5. 1 

Recall that we want to show that Algorithm 5.1 returns a forest with weight that is at least 
a quarter of the weight of the optimal t- restricted forest. Let us distinguish two types of 
constraints: 

(a) the degree of any node is at most t; 

(b) the graph is acyclic. 

Note that the optimal t-restricted forest satisfies both the constraints above, and hence 
the maximum weight set of edges that satisfy both the constraints above has weight at least 
w{F^). Recall that the first stage of Algorithm 5.1 greedily adds edges subject to these two 
constraints — the next two lemmas show that the resulting forest has weight at least ^w{F^). 

Lemma A. 10. The family of subgraphs satisfying the constraints (a) and (b) form a 2- 
independence family. That is, for any subgraph T satisfying (a) and (b), and for any edge 
e E G, there exist at most two edges {ci, 62} in T such that T U {e} — {ei, 62} also satisfies 
constraints (a) and (b). 

Proof. Let T be a subgraph satisfying (a) and (b) and suppose we add e = (m, v) in T. 
Then the degrees of both u and v are at most t + 1. If no cycles were created, then we 
can simply remove an edge in T containing u (if any) and an edge in T containing v (if 
any) to satisfy the degree constraint (a) as well. If adding e created a cycle of the form 
{. . . , (m', u), (u, t>), (t>, t>')}, then the edges {u\ u) and (f , v') can be removed to satisfy both 
constraints (a) and (b). □ 

Lemma A.ll. Let Fi be the forest output after Step 1 of algorithm 5.1. Then w{Fi) > 

\w{f:). 

Proof. Let F** be a maximum weight forest that obeys both constraints (a) and (b). Since 
the optimal t-restructed forest F* obeys both these constraints, we have w{F^) < w{F**). By 
a theorem of Hausmann, Korte and Jenkyns (1980), in a ^-independence family the greedy 
algorithm is a ^-approximation to the maximum weight p-independent set. By Lemma A. 10, 
we know that the set of all subgraphs satisfying constraints (a) and (b) is a 2-independent 
family. Hence, w{Fi) > lw{F**) > |w7(F;). □ 

We can now turn to the proof of Theorem 5.1. 

Proof. Given a graph G, let Fi be the forest output by first step of Algorithm 5.1, and let 
Fa be the forest outputted by the second step. We claim that w{Fa) > |w(Fi); combined 
with Lemma A.ll, this will complete the proof of the theorem. 
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To prove the claim, we first show that given any tree T with edge weights and maximum 
degree t > 2, we can obtain a sub-forest F with total weight w{F) > ^w{T), and where the 
number of edges in each tree in the forest F is at most t — 1. Indeed, root the tree T at an 
arbitrary node of degree-1, and call an edge e odd or even depending on the parity of the 
number of edges in the unique path between e and the root. Note that the set of odd edges 
and the set of even edges partition T into sub-forests composed entirely of stars of maximum 
degree t — 1, and one of these sub-forests contains half the weight of T, which is what we 
wanted to show. 

Applying this procedure to each tree T in the forest Fi, we get the existence of a t — 1- 
restricted subforest F^ C Fi that has weight at least ^w{Fi). Observe that at — 1-restricted 
subforest is a fortiori a fc-restricted subforest, and since w{Fa) is the best t-restricted sub- 
forest of Fi, we have 

w{Fa) > w{F[) > ^w{F^) > ^w{F;), (A.66) 
completing the proof. □ 



A.7.1. An Improved Approximation Algorithm 

We can get an improved approximation algorithm based on a linear programming ap- 
proach. Recall that F** is a maximum weight forest satisfying both (a) and (b). A result of 
Singh and Lau (2007) implies that given any graph G with non-negative edge weights, one 
can find in polynomial time a forest Fsl such that 



w{Fsl) > w{F**) > w{F:), 



(A.67) 



but where the maximum degree in Fsl is t + 1. Now applying the procedure from the proof 
of Theorem 5.1, we get a t-restricted forest whose weight is at least half of w{Fsl)- 
Combining this with (A.67) implies that w{Fsi) > w{F^), and completes the proof of the 
claimed improved approximation algorithm. We remark that the procedure of Singh and Lau 
(2007) to find the forest Fsl is somewhat computationally intensive, since it requires solving 
vertex solutions to large linear programs. 



A. 8. Proof of Theorem 5.2 

Proceeding as in the proof of Theorem 4.2, we have that 

\R{Pf, ) - Rim ) I ^ RiPp, ) - ^-1 (Pf, ) + \Rn, {Pf, ) - R{PFl 

= Op{k(f)nid) + dipnid)) + Rn^{pp^) - R{pf^ 



(A.68) 
(A.69) 



Now, let Hni denote the estimated entropy H{X) = J2k H{Xk), constructed using the kernel 
density estimates Pm (xk). Since the risk is the negative expected log-likelihood, we have using 
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the approximation guarantee that 



RuApfJ - Rim) 



c 

c 

Op (k*Md) + dMd) + ^^w{F:) 



(A.70) 



< 



(A.71) 



(A.72) 



(A.73) 



and the result follows. 

A. 9. The TreePartition Subroutine 

To produce the best t-restricted subforest of the forest Fi, we use a divide-and-conquer forest 
partition algorithm described by (Lukes, 1974), which we now describe in more detail. 

To begin, note that finding an optimal subforest is equivalent to finding a partition of the 
nodes in the forest, where each disjoint tree in the subforest is a cluster in the partition. Since 
a forest contains a disjoint set of trees, it suffices to find the optimal t-restricted partition of 
each of the trees. 

For every subtree T, with root we will find a list of partitions v.P = {v.Pq, v. Pi, ...,v.Pk} 
such that 

1. for i 0, v.Pi is a partition whose cluster containing root v has size i; 

2. v.Pi has the maximum weight among all partitions satisfying the above condition. 

We define v.Pq to be argmax{w(f .Pi), . . . , w{v.Pt)}. The Merge subroutine used in TreePartition 
takes two lists of parititions {v.P,Ui.P}, where v is the parent of Ui, v.P is a partition of 
node V unioned with subtrees of children {ui, . . . , Uj-i}, and Wj.P is a partition of the subtree 
of child Ui] refer to Figure 9. 

Since a partition is a list of clusters of nodes, we denote by Conca.t{v.P2,u.Pk-2) the 
concatenation of clusters of partitions v.P2,u.Pk-2- Note that the concatenation forms a 
partition if V.P2 and u.Pk-2 are respectively partitions of two disjoint sets of vertices. The 
weight of a partition is denoted w{v.P2), that is, the weight of all edges between nodes of 
the same cluster in the partition V.P2. 
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Fig 9: The TreePartition procedure to merge two subproblems. 



Algorithm A. 1 TreePartition(T, t) 

1: Input a tree T, a positive integer t 

2: Returns an optimal partition into trees of size < t. 

3: Initialize v. Pi — [{v}] where v is root of T, if v has no children, return v. Pi 

4: For all children {ui, ...Us} of v, recursively call TreePartition(ui, t) to get a collection of lists of parti- 
tions {ui.P,U2.Pt..Us.P} 

5: For each child Ui G {ui, ...Us} of v 
Update v.P ■(— Merge{ui.P, v.P) 

6: Output v.P[) 



Algorithm A. 2 Merge(t;.P, u.P) 

1: Input a list of partitions v.P and u.P, where is a parent of u. 
2: Returns a single list of partitions v.P' . 
3: For i = 1, . . . , i: 

1. Let [s* ,t*) = argmax(^ u'(Concat(t;.Ps, M.Pt)) 

2. Let v.Pl = CorLC3.t{v.Ps* ,u.Pf) 

4: Select w.Pq = argniax^ p, w{v.Pl) 
5: Output {v.Pi,v.Pi,...v.P^J 



44 



Appendix B: Computation of the Mutual Information Matrix 



In this appendix we explain different methods for computing the mutual information matrix, 
and making the tree estimation more efficient. One way to evaluate the empirical mutual 
information is to use 



Compared with our proposed method 



^ m m ^ I \ 

"^^^ P„i(a;fci)Pni(x£j) 

(B.l) is somewhat easier to calculate. However, if the sample size in D\ is small, the ap- 
proximation error can be large. A different analysis is needed to provide justification of the 
method based on (B.l), which would be more difficult since Pni(-) is dependent on T>x. For 
these reasons we use the method in (B.2). 

Also, note that instead of using the grid based method to evaluate the numerical inte- 
gral, one could use sampling. If we can obtain mi i.i.d. samples from the bivariate density 



(Xr,Xf ) ~ p„,(x„x,), (B.3) 



=1 

then the empirical mutual information can be evaluated as 

^itlt "p(xf))p(A, 



i{x,- X,) = — y log ;/ \\ . (B.4) 



Compared with (B.l), the main advantage of this approach is that the estimate can be 
arbitrarily close to (3.5) for large enough mi and m. Also, the computation can be easier 
compared to Algorithm 3.1. Let p„j(Xj,Xj) be the bivariate kernel density estimator on T>\. 
To sample a point from p„j(Aj, A^), we first random draw a sample ( X]''^ \ xj^ ) from Pi, 
and then sample a point (X, Y) from the bivariate distribution 




(^,^) - T^K ^- K ^- . (B.5) 



Though this sampling strategy is superior to Algorithm 3.1, it requires evaluation of the 
bivariate kernel density estimates on many random points, which is time consuming; the 
grid-based method is preferred. 

In our two-stage procedure, the stage requires calculation of the empirical mutual infor- 
mation I{Xi]Xj) for (2) entries. Each requires 0{rn?ni) work to evaluate the bivariate and 
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univariate kernel density estimates on the mxm grid, in a naive implementation. Therefore, 
the total time to calculate the empirical mutual information matrix M is 0{m?ni <P). In the 
second stage, the time complexity of the Chow-Liu algorithm is dominated by the first step. 
Therefore the total time complexity is O {im?nid?). The first stage requires 0{(P') space to 
store the matrix M and O(m^ni) space to evaluate the kernel density estimates on Vi. The 
space complexity for the Chow-Liu algorithm is 0{(f), and thus the total space complexity 
is 0{(P + m?ni). 

Algorithm B.l More efficient calculation of the mutual information matrix M . 
Initialize M — Odxd and ij'*) = 0„jxm for i = 1, . . . , d. 

% calculate and pre-store the univariate KDE 
for k = 1, . . . , d do 
for fc' = 1, . . . , m do 



vW _ Ji"") 



for fc' = 1, . . . , m do 

% calculate the components used for the bivariate KDE 
for i' ~ 1, . . . ,ni do 
for i = 1, . . . , d do 



h2 \ h2 I 

% calculate the mutual information matrix 
for i' = 1 , . . . , m do 
for i = 1, . . . , d — 1 do 
for j = i + 1, . . . , d do 



for i' = 1, . . . , ni do 

p{xf p{xf \xf ^)/ni 

M(^,J) ^ M(^,J) + J^p(xf\xf) ■ log \ xf ') • 



The quadratic time and space complexity in the number of variables d is acceptable for 
many practical applications but can be prohibitive when the dimension d is large. The main 
bottleneck is to calculate the empirical mutual information matrix M. Due to the utiliza- 
tion of the kernel density estimate, the time complexity is 0{d'^m?ni). The straightforward 
implementation in Algorithm 3.1 is conceptually easy but computationally inefficient, due 
to many redundant operations. For example, in the nested for loop, many components of 
the bivariate and univariate kernel density estimates are repeatedly evaluated. In Algorithm 
B.l, we suggest an alternative method which can significantly reduce such redundancy at 
the price of increased but still affordable space complexity. 

The main technique used in Algorithm B.l is to change the order of the multiple nested 
for loops, combined with some pre-calculation. This algorithm can significantly boost the 
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empirical performance, although the worst case time complexity remains the same. An al- 
ternative suggested by Bach and Jordan (2003) is to approximate the mutual information, 
although this would require further analysis and justification. 
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